function [etrh, detrb, etri, etr2i] = h2gr(params, gin, gtin, gttin)
  Nx = params.Nx;
  Ny = params.Ny;
  Nt = Nx.*Ny;

  bench = true;

  gsteps = size(gin, 1);
  grerr = zeros(gsteps, 1);

  global g0det;

    etrh    = zeros(size(gin, 1), params.Nx, params.Ny);
    detrb   = zeros(size(gin, 1), params.Nx, params.Ny);
    etri    = zeros(size(gin, 1), params.Nx, params.Ny);
    etr2i   = zeros(size(gin, 1), params.Nx, params.Ny);



  % in computing the gr structures we'll loop over the output points, as opposed to compute all of them at the same time using matlab magic in order to converve memory. gamd has 16x times more components!

  % pre-allocate memory
  gd   = zeros(params.Nx, params.Ny, params.Nz, 4, 4, 4); 
  gdt  = zeros(params.Nx, params.Ny, params.Nz, 4, 4, 4); % time derivative of gd
  gamd = zeros(params.Nx, params.Ny, params.Nz, 4, 4, 4, 4);
  Gdd =  zeros(params.Nx, params.Ny, params.Nz, 4, 4);






  for p=1:size(gin, 1)
    % zero out relevant variables
    gd   = zeros(params.Nx, params.Ny, params.Nz, 4, 4, 4);
    gdt  = zeros(params.Nx, params.Ny, params.Nz, 4, 4, 4); % time derivative of gd
    gamd = zeros(params.Nx, params.Ny, params.Nz, 4, 4, 4, 4);
    Gdd =  zeros(params.Nx, params.Ny, params.Nz, 4, 4);

    g  = reshape( gin(p,:,:,:,:,:), params.Nx, params.Ny, params.Nz, 4, 4);
    gt = reshape(gtin(p,:,:,:,:,:), params.Nx, params.Ny, params.Nz, 4, 4);


    % *************************************************************
    % inverse of g, computed as z^2 gi
    % *************************************************************
    [gi git] = h2ginverse(params, g, gt);

    % *************************************************************
    %  t derivative
    % *************************************************************
    % copy time derivatives into gd
    gd (:,:,:,:,:,1) =  gt;
    gdt(:,:,:,:,:,1) = gttin(p,:,:,:,:,:);


    % *************************************************************
    %  x & y derivatives
    % *************************************************************
    for iz = 1:params.Nz
      z = params.Cz(iz);
      for i=1:4
        for j=i:4
          gd(:,:,iz,i,j,2)  = dx(reshape( g(:,:,iz,i,j), params.Nx, params.Ny));
          gd(:,:,iz,i,j,3)  = dy(reshape( g(:,:,iz,i,j), params.Nx, params.Ny));
          gdt(:,:,iz,i,j,2) = dx(reshape(gt(:,:,iz,i,j), params.Nx, params.Ny));
          gdt(:,:,iz,i,j,3) = dy(reshape(gt(:,:,iz,i,j), params.Nx, params.Ny));
        end
      end
    end % end of the z loop
    % fill symmetric
    for i=1:4
      for j=i+1:4
        gd(:,:,:,j,i,2)  = gd(:,:,:,i,j,2);
        gdt(:,:,:,j,i,2) = gdt(:,:,:,i,j,2);
      end
    end


    % *************************************************************
    % z derivative
    % *************************************************************    
    % input g has 1/z^2 factored, and this derivative produces a 1/z^3 
    for i=1:4
      for j=i:4
        gd (:,:,:,i,j,4) = chebgrad(reshape( g(:,:,:,i,j), params.Nx, params.Ny, params.Nz), params.Cz(1), params.Lz, 3);
        gdt(:,:,:,i,j,4) = chebgrad(reshape(gt(:,:,:,i,j), params.Nx, params.Ny, params.Nz), params.Cz(1), params.Lz, 3);
        for (iz = 1:params.Nz)
          gd (:,:,iz,i,j,4) = params.Cz(iz).*gd (:,:,iz,i,j,4);
          gdt(:,:,iz,i,j,4) = params.Cz(iz).*gdt(:,:,iz,i,j,4);
        end
        gd (:,:,:,i,j,4) = gd (:,:,:,i,j,4) + -2.*reshape( g(:,:,:,i,j), params.Nx, params.Ny, params.Nz);
        gdt(:,:,:,i,j,4) = gdt(:,:,:,i,j,4) + -2.*reshape(gt(:,:,:,i,j), params.Nx, params.Ny, params.Nz);
      end
    end
    % fill symmetric
    for i=1:4
      for j=i+1:4
        gd (:,:,:,j,i,4) = gd (:,:,:,i,j,4);
        gdt(:,:,:,j,i,4) = gdt(:,:,:,i,j,4);
      end
    end


    % pull of the overall 1/z^3 factor from the t,x,y derivatives!
    for iz = 1:params.Nz
      z = params.Cz(iz);
      gd(:,:,iz,:,:,1) = gd(:,:,iz,:,:,1).*z;
      gd(:,:,iz,:,:,2) = gd(:,:,iz,:,:,2).*z;
      gd(:,:,iz,:,:,3) = gd(:,:,iz,:,:,3).*z;

      gdt(:,:,iz,:,:,1) = gdt(:,:,iz,:,:,1).*z;
      gdt(:,:,iz,:,:,2) = gdt(:,:,iz,:,:,2).*z;
      gdt(:,:,iz,:,:,3) = gdt(:,:,iz,:,:,3).*z;
    end

    % *************************************************************
    % gd & gdt have now an overall 1/z^3 on all components now
    % *************************************************************    
    


    % *************************************************************
    %  Gamma 1 / z
    % *************************************************************

    % be really careful with the contractions: gi is z^2, g is 1/z^2, gd is 1/z^3
    %   gd - 1/z^3
    %   gi - z^2
    % hence, gam -- 1/z
    for i=1:4
      for j=1:4
        for k=j:4
          gam(:,:,:,i,j,k) = 1/2*( ...
          gi(:,:,:,i,1).*(gd(:,:,:,1,j,k) + gd(:,:,:,1,k,j) - gd(:,:,:,j,k,1)) + ...
          gi(:,:,:,i,2).*(gd(:,:,:,2,j,k) + gd(:,:,:,2,k,j) - gd(:,:,:,j,k,2)) + ...
          gi(:,:,:,i,3).*(gd(:,:,:,3,j,k) + gd(:,:,:,3,k,j) - gd(:,:,:,j,k,3)) + ...
          gi(:,:,:,i,4).*(gd(:,:,:,4,j,k) + gd(:,:,:,4,k,j) - gd(:,:,:,j,k,4)) );
        end
      end
    end
    for j=1:4
      for k=j+1:4
        gam(:,:,:,:,k,j) = gam(:,:,:,:,j,k);
      end
    end




    % *************************************************************
    % Gamma's derivatives
    % *************************************************************

    % time derivative
    for i=1:4
      for j=1:4
        for k=j:4
          gamd(:,:,:,i,j,k,1) = 1/2*( ...
          git(:,:,:,i,1).*(gd(:,:,:,1,j,k) + gd(:,:,:,1,k,j) - gd(:,:,:,j,k,1)) + ...
          git(:,:,:,i,2).*(gd(:,:,:,2,j,k) + gd(:,:,:,2,k,j) - gd(:,:,:,j,k,2)) + ...
          git(:,:,:,i,3).*(gd(:,:,:,3,j,k) + gd(:,:,:,3,k,j) - gd(:,:,:,j,k,3)) + ...
          git(:,:,:,i,4).*(gd(:,:,:,4,j,k) + gd(:,:,:,4,k,j) - gd(:,:,:,j,k,4)) + ...
          gi(:,:,:,i,1).*(gdt(:,:,:,1,j,k) + gdt(:,:,:,1,k,j) - gdt(:,:,:,j,k,1)) + ...
          gi(:,:,:,i,2).*(gdt(:,:,:,2,j,k) + gdt(:,:,:,2,k,j) - gdt(:,:,:,j,k,2)) + ...
          gi(:,:,:,i,3).*(gdt(:,:,:,3,j,k) + gdt(:,:,:,3,k,j) - gdt(:,:,:,j,k,3)) + ...
          gi(:,:,:,i,4).*(gdt(:,:,:,4,j,k) + gdt(:,:,:,4,k,j) - gdt(:,:,:,j,k,4)) );
        end
      end
    end
    for j=1:4
      for k=j+1:4
        gamd(:,:,:,:,k,j,1) = gamd(:,:,:,:,j,k,1);
      end
    end

    % space derivative
    %  for loops only through 1/2 of the symmetric part
    for iz=1:params.Nz
      for i=1:4
        for j=1:4
          for k=j:4 % to compute one part of the symmetric tensor change 1 to j here, and uncomment the portion below
            gamd(:,:,iz,i,j,k,2) = dx(reshape(gam(:,:,iz,i,j,k), params.Nx, params.Ny));
            gamd(:,:,iz,i,j,k,3) = dy(reshape(gam(:,:,iz,i,j,k), params.Nx, params.Ny));
          end
        end
      end
    end
    % fill the symmetric part
    for j=1:4
      for k=j+1:4
        gamd(:,:,:,:,k,j,2) = gamd(:,:,:,:,j,k,2); ;
        gamd(:,:,:,:,k,j,3) = gamd(:,:,:,:,j,k,3); ;
      end
    end

    % pull out an overall 1/z^2 from the t,x,y derivatives:
    for iz=1:params.Nz
      z = params.Cz(iz);
      gamd(:,:,iz,:,:,:,1) = gamd(:,:,iz,:,:,:,1).*z;
      gamd(:,:,iz,:,:,:,2) = gamd(:,:,iz,:,:,:,2).*z;
      gamd(:,:,iz,:,:,:,3) = gamd(:,:,iz,:,:,:,3).*z;
    end

    % compute the r derivative, keeping in mind that gamma has 1/z factored out. Thus gamdz is 1/z^2


    % z derivative 
    for j=1:4
      for k=j:4
        gamd(:,:,:,:,j,k,4) = chebgrad(reshape(gam(:,:,:,:,j,k), params.Nx, params.Ny, params.Nz, 4), params.Cz(1), params.Lz, 3); % the last is dimension

        for (iz = 1:params.Nz)
          gamd(:,:,iz,:,j,k,4) = params.Cz(iz).*gamd(:,:,iz,:,j,k,4);
        end
        gamd(:,:,:,:,j,k,4) = -gam(:,:,:,:,j,k) + gamd(:,:,:,:,j,k,4);
      end
    end
    % fill the symmetric part
    for j=1:4
      for k=j+1:4
        gamd(:,:,:,:,k,j,4) = gamd(:,:,:,:,j,k,4);
      end
    end



    % tic
    % % *************************************************************
    % %  RRRic-ci Scalarrrr
    % % *************************************************************

    % % R = gUU ( Gam Gam + D Gam ) ~ z^2 ( 1/z^2 +1/z^2 ) ~ 1
    % gr = 0;
    % for a=1:4
    %   for b=1:4
    %     for c=1:4
    %       tmp = 0;
    %       for d=1:4
    %         tmp = tmp + gam(:,:,:,d,a,b).*gam(:,:,:,c,c,d) - gam(:,:,:,d,a,c).*gam(:,:,:,c,b,d);
    %       end
    %       gr = gr + gi(:,:,:,a,b).*( gamd(:,:,:,c,a,b,c) - gamd(:,:,:,c,a,c,b) + tmp );
    %     end
    %   end
    % end


    % % *************************************************************
    % %  RUddd
    % % *************************************************************
    % % RUddd ~ D Gam + Gam Gam ~ 1/z^2
    % ruddd = zeros(params.Nx, params.Ny, params.Nz, 4, 4, 4, 4);
    % for a=1:4
    %   for b=1:4
    %     for c=1:4
    %       for d=1:4
    %         tmp = 0;
    %         for k=1:4
    %           tmp = tmp + gam(:,:,:,a,c,k) .* gam(:,:,:,k,d,b) - gam(:,:,:,a,d,k) .* gam(:,:,:,k,c,b);
    %         end
    %         ruddd(:,:,:,a,b,c,d) = (gamd(:,:,:,a,d,b,c) - gamd(:,:,:,a,c,b,d)) + tmp;
    %       end
    %     end
    %   end
    % end


    %  check bianchi identity
    % bianchi = zeros(params.Nx, params.Ny, params.Nz, 4, 4, 4, 4);
    % for a=1:4
    %   for b=1:4
    %     for c=1:4
    %       for d=1:4
    %         bianchi(a, b, c, d) = ruddd(a,b,c,d) + ruddd(a,c,d,b) + ruddd(a,d,b,c);
    %       end
    %     end
    %   end
    % end


    % *************************************************************
    % Rdd
    % *************************************************************
    % 1/z^2

    rdd = zeros(params.Nx, params.Ny, params.Nz, 4, 4);
    for b=1:4
      for d=b:4
        tmp = 0;
        for a=1:4
          for k=1:4
            tmp = tmp + gam(:,:,:,a,a,k) .* gam(:,:,:,k,d,b) - gam(:,:,:,a,d,k) .* gam(:,:,:,k,a,b);
          end
          tmp = tmp + (gamd(:,:,:,a,d,b,a) - gamd(:,:,:,a,a,b,d));
        end
        rdd(:,:,:,b,d) = tmp;
      end
    end
    for i=1:4
      for j=i+1:4
        rdd(:,:,:,j,i)  = rdd(:,:,:,i,j);
      end
    end

    % ********************************
    % einstein tensor
    % ********************************
    edd = rdd + 3.*g;

    % edd -- 1/z^2
    % gi = z^2
    % eUd -- 1

    edU = zeros(params.Nx, params.Ny, params.Nz, 4,4);
    
    for i=1:4
      for j=1:4
        for c=1:4
          edU(:,:,:,i,j) = edU(:,:,:,i,j) + edd(:,:,:,i,c).*gi(:,:,:,c,j);
        end
      end
    end

    % compute the trace:
    etr = zeros(params.Nx, params.Ny, params.Nz);
    for i=1:4
      etr(:,:,:) = etr + edU(:,:,:,i,i);
    end

    % max of abs(tr)
    etrh(p,:,:) = squeeze(max(abs(etr), [], 3));

    % dets
    detedetg = zeros(params.Nx, params.Ny, params.Nz);
    for i=1:params.Nx
      for j=1:params.Ny
        for k=1:params.Nz
          dete = det(squeeze(edd(i,j,k,:,:)));
          detg = det(squeeze(g(i,j,k,:,:)));
          detedetg(i,j,k) = dete./detg;
        end
      end
    end

    detrb(p, :, :) = squeeze(max(detedetg, [], 3));


    

    %% integrate tr root(g) over the manifold
    %etrg = 0.*etr;
    %for zi=1:params.Nz
    %  etrg(:,:,zi) = etr(:,:,zi)./params.Cz(zi).^4;
    %end
    %for x = 1:params.Nx
    %  for y = 1:params.Ny
    %    tmp = chebint(etrg(x,y,:), params.Lz0, params.Lz);
    %    etrh(p,x,y) = tmp(length(tmp));
    %  end
    %end

    % extract the trace at the horizon
    %etrh(p,:,:) = squeeze(etr(:,:,params.Nz));

  
    % compute dtrace
    %detr = chebgrad(etr, params.Cz(1), params.Lz, 3);

    % extract the dtrace at z=0
    %detrb(p,:,:) = squeeze(detr(:,:,params.Nz));

    etr2 = etr.*etr;
    % dummy-integrate
    for x = 1:params.Nx
      for y = 1:params.Ny
        tmp = chebint(etr(x,y,:), params.Lz0, params.Lz);
        etri(p,x,y) = tmp(length(tmp));

        tmp = chebint(etr2(x,y,:), params.Lz0, params.Lz);
        etr2i(p,x,y) = tmp(length(tmp));
      end
    end



    refresh
    pause(0.001)

  end % end of master loop over output t points



 

end
